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ABSTRACT 

Using the largest cosmological reionization simulation to-date {^24 billion particles) , we use the 
genus curve to quantify the topology of the neutral hydrogen distribution on scales > 1 Mpc as it 
evolves during cosmological reionization. We find that the reionization process proceeds primarily in 
an inside-out fashion, where higher density regions become ionized earlier than lower density regions. 
There are four distinct topological phases: (1) Pre- reionization at z > 15, when the genus curve is 
consistent with a Gaussian density distribution. (2) Pre-overlap at 10 < z < 15, during which the 
number of HII bubbles increases gradually with time, until percolation of HII bubbles starts to take 
effect, characterized by a very flat genus curve at high volume fractions. (3) Overlap at 8 < z < 10, 
when large HII bubbles rapidly merge, manifested by a precipitous drop in the amplitude of the genus 
curve. (4) Post-overlap at 6 < z < 8, when HII bubbles have mostly overlapped and the genus curve 
is consistent with a diminishing number of isolated neutral islands. After the end of reionization 
(2 < 6), the genus of neutral hydrogen is consistent with Gaussian random phase, in agreement with 
observations. 

Subject headings: cosmology: theory — intergalactic medium — large-scale structure of universe — 
methods: numerical — methods: statistical — radiative transfer 



1. INTRODUCTION 

Two recent major observational milestones combine to 
suggest that the so-called "dark ages" of the universe 
(jPeacock 1992t iRees 19981) may be rather complex. The 
first major milestone was laid down by the recent absorp- 
tion spectrum observations of high redshift quasars from 
the Sloan Digital Sky Survey (SDSS), indicating that the 
final reio nization episode came to completion a t z ~ 6 
(seei; i; iBarkana 200l iCen fc McDonald 20021 : IH). 

The Wilkinson Microwave Anisotropy Probe (WMAP) 
polarization observations, probing the ionization state of 
intergalactic medium (IGM) in the high redshift universe, 
set the second major milestone. The WMAP observa- 
tions indicate Te = 0.09 ± 0.03 ()23l ). suggesting that the 
process of reionization occured at redshifts considerably 
greater than z ~ 6, possibly z > 15, although the current 
uncertainty in Tg allows a very wide range of reionization 
histories. 

Although the complex reionization process implied by 
these observation s did not appear as a to tal surprise (e.g. , 
Gncdin 2000b; IBarkana fc Loeb 2001] : jMadau 2002 : 
Wvithe & Locb 2003": '16': 'Cen 2003; '25';_[22 g 
Venkat esan & Truran 2003; Ricotti & Ostri ker 2004( 1. it 
is unclear at this time how the universe had actually 
evolved from z ~ 1000 to z = 6, for different scenarios 
could be designed to arrive at the same observed total 
Thomson optical depth (e.g., |l3) as well as satisfy the 
sudden rise in ionizing radiation background observed at 
z~ 6. 

Thus, at present, many details of the cosmological 
reionization process are still unknown. It is expected, 
however, that upcoming observations at various obser- 
vational bands from radio to optical promise to provide 
much more refined statistical information than what is 
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available today. Therefore, advances on the theoretical 
side are demanded in order to provide both useful in- 
sights and a quantitative framework for the interpreta- 
tion of future observations. 

At present, we make the first attempt to quantify the 
topological evolution of cosmological reionization, in the 
context of the conventional scenario that stellar radia- 
tion is primarily responsible for reionizing the universe. 
We utilize a state-of-the-art cosmological reionization 
simulation of a box size 100 Mpc which involves a 
hybrid algorithm i nclusive of detailed radiative transfer 
()Trac &: Cen 2006t [2l[ ). This paper is organized as fol- 
lows: the method section describes the genus and simu- 
lation code, followed by a description of the results and 
the manner in which they are presented. We then give a 
detailed discussion of the changes in the topology of the 
neutral fraction at different stages of reionization. 

2. METHOD 
2.1. Genus: A quantitative measure of topology 

In order to study the topology of a three-dimensional 
data set like those in the simulation, we must first spec- 
ify two-dimensional contours. From a continuous den- 
sity distribution, this can be done by defining contours 
of constant density at different volume fractions /, i.e. 
the fraction of the volume occupying the high-density 
side of the contour. The topology of contours at a given 
threshold density can then be quantified with the genus. 

For a Gaussian distribution. 
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(jjsl), where v = 6/^(0)^^'^ is the number of standard de- 
viations by which the threshold density S of the contour 
departs from the mean density in a Gaussian distribution 
with covariance function ^(r). For example, the density 
contour with = 1.5 is that which encloses ~ 93% of the 
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volume, or f — 0.07 (i.e. 7% of the volume has a higher 
density than v — 1.5). 

While the relation between / and v as defined in Equa- 
tion [1] holds only for a Gaussian distribution, it has be- 
come the convention to use z/ as a convenient label for 
the different density contours when studying the genus 
of a distribution, even if the distribution is not Gaus- 
sian. In order to allow comparisons with the published 
literature dealing with the genus, we will mostly follow 
this convention of referring to contours using as a label, 
which is related to the actual volume fraction / by Equa- 
tion [H This also has the advantage that large values of 
V intuitively correspond to high threshold densities and 
vice- versa, whereas / decreases with increasing threshold 
density. 

Mathematically, the topology of an object is de- 
fined by its homeomorphism, i.e. its geometric proper- 
ties under deformations and not transformations which 
'brea k' an object or 'conne ct' it with others (see, 
e.g. ISeifert fc Threlfall 1980[ ). One quantitative mea- 
sure of an objects topology is its genus (^13; 8; 11). We 
define the genus of a repeating or infinite contour as 

Gs = (no. of holes) — (no. of isolated regions) (2) 

where a 'hole' is like that in a torus, while an 'isolated 
region' can be either above or below the threshold den- 
sity of a contour In this paper we define 'holes' 
and 'isolated regions' depending on whether the thresh- 
old density of the contour being studied is higher or lower 
than the 50% median density value: for voids, isolated 
regions are defined as those with lower densities than 
the threshold density, while the reverse is true for dense 
pockets. 

Thus, a single dense region in an empty box and a 
single void in a solid box are both counted as 'isolated 
regions' and are hence topologically equivalent. This al- 
lows for a symmetry in the definition of the genus for low 
density [v < 0) and high density {i/ > 0) regions. 

The genus of a repeating or infinite two-dimensional 
contour can be calculated by integrating its Gaussian 
curvature K over the entire area A, 

= ^ / ™ (3) 

where the Gaussian curvature is defined as the reciprocal 
of the two principle radii of curvature oi and 02 at the 
point, i.e. K = 1/(0102) (jHHH). As an illustration, take 
for example a density contour which comprises of 20 iso- 
lated spheres, each of radius r^. The Gaussian curvature 
at any point on the surface of each sphere is K = rj^, 
while the total surface area is of course 20(47rr^). Using 
Equation [3] gives a genus of -20 as we would expect from 
the definition given in Equation[21 A multiply-connected, 
sponge-like contour would have a positive genus value, 
since the two principal radii of curvature on its surface 
tend to point in opposite directions rather than in the 
same direction, giving it a negative curvature. In a three- 
dimensional data set where each data pixel is essentially 
a cube, the genus is calculated by summing up the angle 
deficits at the vertices of the Lego-like contour surface. 
We use a c omputer program , ContourSD, written by 
(author?) (| Weinbe rg 19881 for this purpose. 

A point worth noting is that since the genus is calcu- 
lated over the entire volume of the box, a combination of 



multiply-connected structures and isolated pockets could 
lower the amplitude of the genus. We thus have to be 
careful in distinguishing between, say, a volume with one 
figure-8 shaped object plus three small isolated spheres 
in it, as opposed to another volume with just one large 
isolated sphere, since both situations will give a total 
genus of -1. 

The genus provides a useful tool for probing the rela- 
tive connectedness of high- and low-density regions in the 
large-scale structure of the universe, for scalar quantities 
such as the dark matter distribution or density of HI. In 
addition to being dependent on the underlying distribu- 
tion, the topology of large-scale structure also depends 
on the density threshold specified for a given contour sur- 
face. In general, contours of low threshold densities (i.e. 
voids occupying lowhigh volume fractions) have a differ- 
ent topology from contours of high threshold density (i.e. 
dense regions occupying low volume fractions) within a 
given data set. This provides an additional degree of free- 
dom with which to study the topology of a density distri- 
bution, since we can now study the systematic variation 
of Gs{v) with threshold densities using 'genus curves'. 

In addition, the genus is an useful cosmological statistic 
for another reason. The standard big bang-inflationary 
model predicts that small-amplitude density fluctuations 
arise from a random phase Gaussian distribution (2; 1). 
For such a distribution, the mean genus per unit vol- 
ume Qs ^ Gs/V is given by a simple analytic expression 
(IDoroshkevich 19701 : lAdler 19811 : lltl: 

gs=N{l-v^)e--"'^ (4) 

where the normalization N as function of wavenumber k 
is given by 

47r2 V 3 ) 47r2 \ P'{k)d\ ) ^ ' 

where P'{k) — P{k) e^"^ /^''sm is the power spectrum 
which has been smoothed with a smoothing radius rgm- 
We thus see that any Gaussian random-phase density dis- 
tribution would have a curve of gsii^) with the same dis- 
tinctive shape regardless of the power spectrum, which 
enters only through the amplitude. The implication is 
that the shape of the genus curve would remain un- 
changed through linear evolution of the fluctuations. The 
initial conditions of our simulation at z 25 display this 
characteristic curve (top left panel. Figure [2]), and it is 
virtually indistinguishable from the theoretical curve . 

Equation [4] also shows that there is a symmetry be- 
tween positive and negative v in a random phase distri- 
bution, i.e. high- and low-density regions with the same 
absolute value of v would have an identical topology. For 
\u\ > 1 (i.e. / < 16% or / > 84%) , the high- and low- 
density regions would comprise of isolated dense pockets 
and voids respectively as indicated by the negative genus. 
The median density contour with v — or f — 0.5, on 
the other hand, has a positive genus indicating a sponge- 
like, multiply connected structure (f^: [III). 

Historically, the genus has provided strong evidence in 
support of the Gaussian random phase hypothesis (and 
hence the now-standard inflationary big-bang model). 
Studies on the topology of large-scale structure as ob- 
served in galaxy surveys have yielded genus curves in 
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remarkable agreement with the form shown in Equation 
m apart from small-scale deviations due to non-linear 
grav itational effects and biasing from galaxy formation 
d; Ei II I)- In the present context, this means 
that the genus is a powerful and sensitive probe of the 
topology of the neutral fraction in the universe, as reion- 
ization causes a departure from a Gaussian random phase 
distribution. 

2.2. Simulations 

We use th e largest simulatio n of cosmic reionization 
run to date (|Trac fc Cen 2006t EU). It was run with a 
hybrid code containing a N-body algorithm for dark mat- 
ter, prescriptions for baryons and star formation, and a 
radiative transfer (RT) algorithm for ionizing photons. 
Here we summarize the main simulation parameters. 

The hybrid simulation was run with the cosmologi- 
cal parameters: f2„ = 0.26, Qi = 0.74, = 0.044, 
h = 0.72, as = 0.77, and Ug — 0.95, based on the lat- 
est resuhs from WMAP, SDSS, BAG, SN, and HST (see 
? , and references therein). In a lOOMpcft."^ simula- 
tion box, a high resolution N-body calculation for 2880'^ 
(i.e. 23.9 billion) dark matter particles on an effective 
grid with 1 1520^ cells is perfor med using a particle- multi- 
mesh code (jTrac fc Pen 20061 ) . With a particle mass res- 
olution of 3.02 X 10^ Mq/i-i, halos can be reliably re- 
solved down to masses of ~ IQ^MqH^^ , accounting for 
the majority of photoionizing sources. 

The radiative transfer of ionizing radiation was run si- 
multaneously with the N-body calculations using a RT 
grid with 360"^ cells. However, the ionization and recom- 
bination calculations were done for each particle indi- 
vidually rather than on the grid to preserve small-scale 
information down to scales of several comoving kpc h~^. 
For post-processing, the dark matter, baryons, and ra- 
diation are collected on a grid with 720^ cells and the 
data is saved every 10 million years from z — 25 down to 
z = 5. 

3. RESULTS 

In the present work, we study the density distribu- 
tion of neutral hydrogen (HI) at various redshifts dur- 
ing the epoch of reionization. The redshift evolution of 
the volume- and mass-weighted HI fractions are shown 
in Figure [1] In our simulation, we can not resolve the 
damped Lyman alpha systems associated with cooled 
neutral gas in halos above the cooling mass. Further- 
more, we can only account for Lyman limit systems in the 
form of neutral gas in partially self-shielded mini-halos 
down to about half a dex in mass below the cooling mass. 
Thus our HI densities are representative of the IGM and 
partially exclusive of collapsed regions. We have chosen 
values for the star formation efficiency and radiation es- 
cape fraction such that complete reionization is achieved 
by z ~ 6, as suggested by observations of Lyman alpha 
absorption in high redshift quasars. The results of our 
analysis can also be applied to other reionization histories 
since recent work ()28l : IITI ) have shown that the distribu- 
tion of HI is, in general, similar for different reionization 
histories when compared at the same ionization fraction. 
We thus use the volume-weighted ionization fraction /hi 
in our plots as an alternative label for different epochs in 
addition to the redshift z. 

The lOOMpc/i^^ simulation box has 720^^ cells, i.e. a 
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Fig. 1. — The volume- weighted (soUd) and mass- weighted (dot- 
ted) mean HI fractions from our simulation plotted as a function 
of redshift. These enable our results to be applied to other reion- 
ization histories by matching the mean HI fractions at different 
epochs. 

spatial resolution of ^ 0.14Mpc/i~^. Before beginning 
the analysis, we first smooth the data via a Fourier trans- 
form with a Gaussian kernel W = /^''sm. We use a 
smoothing length of rgm — lMpc/i~^, a scale which is 
adequate to smooth over individual pixels without wash- 
ing out the structure of HII bubbles during reionization. 

We apply the Contour3D code to the simulation 
results at redshifts ranging from z '--^ 25, well before 
the reionization process, to z ~ 5.5, after the end of 
reionization. Figure O shows the genus curves of the HI 
distribution at 9 different redshifts during reionization. 
These particular redshifts are shown as they represent 
the 'turning points' in the evolution of the topology, and 
the change between any two adjacent redshifts plotted in 
Figure [2] is approximately monotonic. 

Before we present the detailed discussion, we note that 
the actual HI threshold density pthres,Hi corresponding to 
a given volume fraction will not, in general, be constant 
as the universe evolves. As there is a spatial dependence 
on when neutral regions get reionized, an initially under- 
dense HI region could very well become a dense region 
relative to the rest of the reionized universe at a later 
time. We thus need to keep track of the redshift and 
spatial variation of the threshold density values in order 
to correctly interpret the changes in topology. 

For example, in the initial conditions at z = 25 be- 
fore the onset of reionization, the median HI density at 
i/ = 0or/ = 0.5is essentially the same as the overall 
mean hydrogen density , i.e. pthros, m/ipn) ~ 1 (re- 
call that for a Gaussian distribution, v is the number 
of standard deviations from the mean, so the median 
50% volume fraction threshold density is the same as the 
mean density at v = 0) . By the end of reionization 
at z ^ 6, the average neutral fraction is close to zero, so 
Pthros, hi/ (pb) ~ 0. While we expect the densities in gen- 
eral to fall as reionization progresses, from Figure [3] it is 
obvious that /9thrcs,Hi for different values of v decrease at 
different rates during reionization, and at different red- 
shifts. We thus need to keep track of the evolution and 
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Fig. 2. — Genus curves of the HI density distribution at various epochs during reionization, labeled by redshift z and corresponding 
volume- weighted mean neutral hydrogen fraction /hi- Following convention, we label the bottom axes with the quantity v which is related 
to the volume fraction / (on top axes) by Equation [T] (see discussion in text). The dotted curves plotted in the top-left and bottom-right 
panels are the best-fit theoretical curves for a Gaussian random-phase distribution (Equation |4]l . The vertical and horizontal dashed lines 
provide fiducial points for u = and Gs = respectively. 



spatial variation of the actual threshold density values 
corresponding to a given ly in order to correctly interpret 
the changes in topology. 

Another property we can study with the genus analysis 
is the characteristic scale of the neutral structures at a 
given epoch. From Equations H] and [51 the amplitude 
iV of a Gaussian distribution's genus curve is related to 
the cube of its root-mean-squared (r.m.s.) wavenumber 
{k^y^^. We thus have an estimate of the characteristic 
scale of a Gaussian distribution by inverting Equation [5] 
to calculate (A:^)^/^ from the amplitude Gs,poak of the 
genus curve. 

Even if the distribution is not exactly Gaussian 
random-phase, the maximum amplitude of the genus 
curve can still give information about its characteristic 
scale at higher densities so long as the distribution re- 
mains approximately Gaussian around those higher den- 
sities. As an analogy to illustrate this argument, consider 



a sea sponge immersed in water. When the sponge is 
moved to dry land, the topology of the enclosed volume, 
evaluated at the density of water, will change drastically 
as water drains and evaporates away. At the density of 
the sponge itself however, there will be little change in 
the topology compared with when it was immersed in 
the sea. Similarly, even as ionized bubbles begin grow- 
ing in the simulation box, we can still use the amplitude 
of the higher density regions as an indication of scale so 
long as the bubbles do not occupy a major fraction of 
the volume. 

A cursory glance at Figure [2] shows that even when the 
topology deviates from the W-shaped Gaussian curve, if 
we ignore the negative v (i.e. low-density) part of the 
genus curves, the high-density part does indeed retain 
an approximately Gaussian shape down to ^ ~ 10. For 
example, the i' ^ —I portion of the z — 12.1 curve looks 
like part of an off-center Gaussian curve. Immediately 
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Fig. 3. — (a)Left panel shows the contour threshold densities as a function of u, (or equivalently volume fraction /; see Equation [TJ , at 
different redshifts. Note that the threshold densities do not decrease monotonically in the early phases, between 2 ~ 25 and 2 ~ 17 (see 
discussion in text) (b) Right panel shows the redshift evolution of the threshold densities at / = 0.93 or = —1.5 (dotted line), / = 0.5 or 
u = Q (solid line) and / = 0.07 or v = 1.5 (dashed line). 

4.1. Pre-reionization; z > 16; /hi 0.96 

As the pre-ionization universe is completely neutral at 
z ^ 25, the topology of the HI distribution is that of the 
primordial Gaussian random phase distribution (top-left 
panel of Figure [2]) • The HI density contours of this dis- 
tribution are illustrated in the top-left panel of Figure [D 
(for clarity, we only show a 25 Mpc hr^ section extracted 
from the overall 100 Mpc box). Both high- and low- 
density contours are split into similar numbers of isolated 
regions corresponding to a negative genus, while the me- 
dian density contour at i/ = or / = 0.5 is essentially one 
multiply-connected object with many holes in it, giving 
a positive genus. 

While star-formation is already occurring dur- 
ing this epoch, the star-formation rate SFR < 
10^"^ Mq Mpc^'^yr^^ is too small to provide significant 
ionizing flux to the IGM; we see from Figure [1] that 
jm ~ 1- Hence, the dominant effect on the neutral gas 
distribution in this era is linear gravitational evolution. 

In Figure [3l we see that throughout this epoch the 
HI threshold density ati^ = 0or/ = 0.5 remains al- 
most the same as the mean density of all hydrogen in 
the universe, pthres, hi/(ph) ~ 1- Indeed, as the uni- 
verse evolves during this time, the densities of regions 
with v > 1 actually increase relative to the early con- 
ditions at z 25 while the densities of the voids with 
v < 1 decrease slightly. Both curves however remain 
straight lines. This behaviour is as expected from linear 
evolution: the dense regions cluster due to gravitational 
interactions at the expense of the underdense voids. This 
leaves the topology essentially unchanged, as can be as- 
certained by comparing the plot for z — 17.2 in Figure[2] 
with that of the original conditions at z = 25: they are 
almost identical. 
The r.m.s wavenumber of the HI structures also re- 




FlG. 4. — The r.m.s wavenumber of the HI distribution, plotted 
as a function of redshift. The right axis labels f the corresponding 
genus amplitudes Gg p^ak foi' our simulation volume (see Equation 
O. We truncate the plot at around the epoch when the topol- 
ogy departs strongly from a Gaussian distribution, at 2 ^ 8 (see 
discussion in text). 

after z ~ 10, the topology of the distribution becomes 
completely non-Gaussian (e.g. z = 8.2 curve in Figured]) 
and it becomes difficult to define the amplitude of the 
peak. We plot the genus amphtude Gs.pcak — NV (for a 
box volume V = 100'^ Mpc'^/i^^) and the corresponding 
(fc2)i/2 at redshifts 9 < z < 25 in Figure H 

4. DISCUSSION 

In this section, we present a detailed analysis of the 
different phases of reionization identified from the varia- 
tion in the topology of the HI density. In particular, we 
have identified processes corresponding to 'pre-overlap', 
'overlap' and ' post-overlap' r eionization first described 
by (author?) (|Gnedin 2000af) . 



mains approximately constant at {k 



2\l/2 



0.7 Mpc" 



(corresponding to a characteristic scale (A^)^/^ ~ 
9Mpc/i^^) during this epoch as can be seen in Figure 
131 Again, this is characteristic of linear evolution, which 
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<f [HI)> = 1.0 




<f (HI)> 




= S.17 

<f [HI)> = 0.2 9 




= 5.57 
<f(HI)> = 0.0008 



Fig. 5. — HI density contours from a 25Mpc/i~^ subsection of the overall 100Mpc/i~^ simulation box, labeled by redshifts and their 
corresponding volume-weighted mean neutral fractions. The different colors correspond to densities occupying different volume fractions 
at that redshift. Blue: / = 0.93 or = —1.5; green: / = 0.5 or u = 0; red: / = 0.07 or v = 1.5. Note that actual threshold density values 
vary with redshift (see discussion in text). 



does not change the topology in any way. 

At the end of this phase, we begin to see a rise in the 
amphtude of the genus (see plot for z — 15.7 in Figure^]). 
This is an indication that the reionization is becoming 
more significant as HII bubbles begin to grow, leading to 
the next phase. 

4.2. Pre-overlap; 10 < z < 16; 0.63 < /hi < 0.96 

Starting at a redshift of z ~ 15, the topology starts to 
change appreciably as reionization begins in earnest. 

The genus curves (e.g. the plot for z = 15.7 and z = 
12.1 in Figure[2]) start to display an asymmetry between 
low- and high-^ portions: while the absolute value of 
the genus for intermediate and high density regions (i^ > 
0) increase along with the overall amplitude, the curve 
flattens out at low i/, and the absolute genus values at 
these densities decrease. 

A genus curve with this shape is not new: when Gott 
(ill ) studied the topologies of non-random phase distri- 



butions for the universe, one possibility they considered 
was one in which galaxies are distributed at surfaces sur- 
rounding large voids. This has striking similarities with 
the pre-overlap genus curves in the present case, indicat- 
ing clearly that large ionized HII bubbles are being cre- 
ated during this epoch. However, in the present work, 
the bubbles have continuously varying HI densities in 
contrast to the uniform voids studied in Gott (fill). In- 
deed, the increase of the absolute value of the genus to- 
wards the median [v = 0) in the v < regime is in- 
dicative of the intuitive fact that the most highly ionized 
bubbles are rarer than less ionized ones, thus less bub- 
bles are detected at the lowest threshold densities. Also, 
the topology at high densities is different in both cases. 
In Gott (11), the galaxies were distributed randomly on 
the bubble walls (and hence approximately uniformly), 
while in the present case of pre-overlap reionization the 
high density regions outside of the ionized bubbles still 
have an approximately Gaussian distribution of neutral 
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hydrogen. 

From Figure we see that the curves relating thresh- 
old density pthros, hi/ (ph) to a given volume fraction / or 

V have a different shape during this epoch, with a straight 
line at 1/ > 0, but with a tail dropping off steeply towards 
low HI densities at low a characteristic which becomes 
increasingly marked with decreasing redshift. This can 
also be seen in Figure [3|3, in which the threshold den- 
sity for low volume fractions begin decreasing rapidly at 
around z 14, while the threshold densities at median 
and high volume fractions decrease more gently. This 
is indicative of the drastic reduction in neutral density 
within the interiors of HII bubbles. 

In Figure [Sj a comparison between the contours at 
z = 24.9 and z — 9.6 reveals that regions initially oc- 
cupied by high density contours have now become low 
density voids. For example, it is clear that the large high- 
density pocket near the middle front of the box in the 
initial conditions has become a low-density void. This 
is a strong indication of 'inside-out' reionization. In this 
picture, the highest density regions are the first in which 
halos can accrete sufficient mass (and hence virial tem- 
peratures of T > lO**) for atomic cooling to occur and 
cause the halos to collapse and commence star forma- 
tion. These star-forming regions in turn becomes sources 
of ionizing flux, which ionize the surrounding neutral gas. 

This effect means that the HII bubbles detected in the 

V < part of the genus curve are in high-density regions 
which have been strongly ionized, and not regions with 
intrinsically low density. 

During this stage, the ionizing fluxes arc shielded from 
the IGM by virtue of being within dense regions with 
high recombination rates, hence the universe outside of 
the HII regions is little affected by reionization. Outside 
of the increasingly low HI densities in the HII bubbles, 
the relative densities are still the same. For example, a 
ly = —1.5 void in the z = 24.9 distribution is by defini- 
tion less dense than a = median density region during 
the same era, but they still have approximately the same 
relative densities with respect to each other during the 
pre-overlap stage of reionization as they experience lit- 
tle ionizing flux. The main difference is that the original 
dense pockets have become HII regions and are now voids 
with respect to the rest of the universe, so the volume 
fractions covered by the un-ionized regions have been re- 
duced (i.e. v increased). This can be seen by looking 
at the positions of the intermediate (green, = 0) and 
dense (red, v = 1.5) density contours for z = 9.6 in Fig- 
ure [5l the former are in regions formerly occupied by 
the voids while the latter are in regions once part of the 
sponge-like median density contour. The correspondence 
between the contours in is not exact because the actual 
threshold density values for e.g. the ly = 1.5 contour at 
z — 9.6 and ly — a,t z — 24.9 are not exactly matched 
to each other, but is sufficient to demonstrate the broad 
trends at different densities. 

This argument is supported by the evolution of the 
threshold densities with redshift in Figure [3^. Wc sec, in 
this context, that the low-density tail signifying HII bub- 
bles must be from regions that used to be in the high-i^ 
part of the curve in earlier times. The straight portion 
at > is from regions that once occupied the low den- 
sity ^ portion of the straight line representing the 
Gaussian distribution before reionization began. These 



have now been translated to occupy higher volume frac- 
tions relative to the overall distribution. However, even 
in these regions, some small amount of ionization is tak- 
ing place at scales below the smoothing length, which 
slowly decrease the local HI densities. 

As HII regions begin forming in the regions of highest 
baryon density, recombination initially mitigates the ef- 
fects of the ionizing flux. Thus, these partially ionized 
regions initially appear as bubbles of intermediate neu- 
tral density growing within highly dense regions. This 
has the effect of creating holes and isolating individual 
dense clumps within the dense regions as the HII regions 
appear, breaking up the topology and increasing the am- 
plitude of the genus. 

As reionization progresses, the absolute genus value of 
the low-density regions decreases due to the merger of 
small HII bubbles into a smaller number of large bubbles 
(compare v < —1 tails of z = 12.10 and z = 9.6 in 
Figure [2|). This is in agreement with Shin (pll). which 
found that the majority of bubble mergers occur during 
this epoch. 

4.3. Overlap; 8<z< 10; 0.26 < /m < 0.63 

At z '--^ 9, volume- weighted mean HI fraction of the uni- 
verse drops to /hi ^0.5 (Figure [1]), in contrast with the 
leisurely decrease in previous phases, even though there 
is merely a modest increase in the amount of ionizing flux 
emitted by the star-forming r egions (now domin ated by 
Pop II stars; see (author?) (|Trac fc Cen 20061 ) ). This 
sudden decrease occurs as the ionized bubbles overlap 
with each other, so that ionizing photons can now stream 
between previously isolated bubbles and increase the ion- 
izing flux experienced by points within the overlapping 
ionized regions. This increased flux overcomes the effects 
of recombination and increases the ionization rate. At re- 
gions occupying the lowest volume fractions, the gas is es- 
sentially completely ionized, with pthros, m/ipn) ^ 0.01 
at = -1.5 (Figure [3b). 

As the regions occupying the median volume frac- 
tions (/ — 0.5 OT i> ^ 0) are now also part of signifl- 
cantly ionized regions, the genus of what were previously 
multiply-connected regions is now strongly suppressed 
due to the overlapping of HII regions, which closes gaps 
in the sponge-like contours (and hence reduces the pos- 
itive genus). By z = 9.6, the overall amplitude of the 
genus curve is considerably lower than that at z = 12.1 
(Figure [2]). This reflects an increase in the characteristic 
scale of the neutral distribution, driving down the r.m.s. 
wavenumber and hence the genus amplitude (Figure 2]). 
Immediately after z ~ 9, the topology of the distribu- 
tion across all densities becomes strongly non-Gaussian 
as can be seen from a comparison of the genus curves 
for z = 9.6 and z = 8.2 in Figure [D After this point, 
there no longer exists large neutral structures with ap- 
proximately Gaussian distributions within them. This 
means we can no longer estimate a characteristic scale 
using Equation [5] and so we truncate Figure |4] at z ~ 9. 

The overlap of ionized regions means that the topology 
at low-densities is dominated by a small number of large 
interconnected regions, which is detected as > in 
the 1/ < portion of the genus curve for z = 8.2 (Figure 
[2]), and is seen in the corresponding contours in Figure 

El 

With the percolation of HII bubbles through the IGM, 
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the remaining regions with a relatively high HI density 
(yOthres, Hi/ (ph) ^ 0.5 at > 0) are now isolated into indi- 
vidual pockets (see the corresponding contours in Figure 
[51). This is also seen in the high-density j/ > 1.5 portion 
in the genus curve for z — 8.2 in Figured which retains 
vestiges of its curved Gaussian shape. This is because 
these high density regions are still shielded from the ion- 
izing flux flooding the rest of the universe, thus within 
them the relative density distributions still reflect those 
of the initial conditions. But even within these pockets, 
small amounts of star formation are going on at scales be- 
low the smoothing length, which continue pushing down 
slowly on the densities. The stage is now set for the 
reionization of the remaining large-scale HI structures in 
the universe. 

4.4. Post-overlap; 6 < z < 8; /hi < 0.26 

By redshifts of z ^ 8, the median HI density of universe 
is Pthres, Hi/(/3H> ^ 0.1 and dropping (Figure [p. The 
contours with threshold densities below the median, v < 
0, are now almost thoroughly ionized, as characterised 
by the flat tails in the corresponding pthrcs, hi/(ph) v.s. 
ly plots in Figure [3^. 

In the genus curve for z = 7.3 (Figure [2]), the genus at 
V < has now assumed a shape that appears similar to 
that of a Gaussian distribution, but truncated where it 
crosses the median volume fraction at « 0. This in- 
dicates that regions below these threshold densities are 
totally ionized, with the unresolved neutral regions (i.e. 
galaxies) within them acting as tracers of the baryon's 
Gaussian random phase density distribution. A compar- 
ison with the curve in Figure [3] at this redshift concurs 
with the notion that a straight line in the pthres, hi/ (ph) 
v.s. v curve is indicative of a random phase distribution. 
At < 1/ < 1.5, ionized regions have yet to be com- 
pletely ionized, and have the topology of a unified region 
punctuated by small numbers of holes. 

These holes are caused by the remaining resolved pock- 
ets of high-density HI in the universe, which by this epoch 
has a density about 2 orders of magnitude greater than 
fully ionized regions occupying most of the IGM (Fig- 
ure [3)3) . These pockets can also be seen in the genus 
curve at 1/ > 1.5. However, the small negative amplitude 
of the genus indicates that only a small number of neu- 
tral pockets remain, and corresponding volume fraction 
they occupy is small (/ — 0.02 at 1/ — 2). 

With the universe now largely ionized, most of the IGM 
has a low optical depth to ionizing photons. The re- 
maining HI regions are isolated pockets occupying small 
volumes, which implies that they have relatively large 
surface areas exposed to large amounts of ionizing flux. 
Under these circumstances, the final reionization of these 
pockets happen within a relatively short period of time as 
illustrated in Figure [3^, where pthrcs, hi/(ph) at ~ 1.5 
decreases by about 2 orders of magnitude during the in- 
terval 5.6 ^ z <7. 

The genus curve at 2; = 6.4 (Figure[2]) offers an interest- 
ing snapshot of the topology during this time. The curve 
at 1/ < 1 is distinctively that of the Gaussian random- 
phase distribution, corresponding to the baryonic distri- 
bution traced by the unresolved neutral HI. At 1/ > 1, 
we see that the densest ^ 10% (/ ~ 0.1) of the volume 
are in a very tiny number of dense neutral pockets, but 
they no longer exist by z = 5.6. 



At the end of post-overlap reionization at z = 5.6, the 
HI threshold density is a straight line with respect to v 
(Figure [3^), with pthrcs, hi/(ph) ^ 0.01. This indicates 
that the reionization of the universe is complete by this 
redshift, and the tiny amounts of neutral HI remaining 
in the universe now traces the overall baryon distribution 
in the universe, which is still Gaussian random phase to 
a good approximation. This is supported by the corre- 
sponding genus curves (Figure [2|) and contours (Figure 

We note that a post-reionization Gaussian random- 
phase distribution is directly supported by observational 
evidence: Weinberg ([27l ) had studied the 1-dimensional 
topology of neutral hydrogen using Ly-a absorption lines 
in quasars with z < 2.5, and it was found that the dis- 
tribution of threshold crossings as a function of (i.e. 
number of threshold crossings oc exp— 1/^/2) smoothed 
on scales of lMpch~^ is consistent with a Gaussian 
random-phase distribution . 

In the contours plots, the distribution of the high- and 
low-density regions at z = 5.6 bear a resemblance to 
the initial conditions at z = 24.9. For example, the 
anvil-shaped void at the top vertex facing the reader can 
clearly be matched with an similarly-shaped object in 
the initial conditions, as is the small void at the lower 
vertex. In general, the contour distributions in the ini- 
tial and post-reionization volume can be approximately 
matched to each other, but there are small differences. 

The post-reionization universe is in general choppier 
than the initial conditions, most noticeably in the dense 
regions. This is supported by a comparison between the 
two corresponding genus curves (Figure [2]) . The ampli- 
tude of the genus curve at z — 5.6 is about twice that 
at z = 24.9, giving a corresponding decrease in the char- 
acteristic wavelength according to Equation [S] This dif- 
ference is most likely due to non-linear effects which oc- 
curred in the baryon distribution whilst reionization was 
taking place, but is beyond the scope of this paper. 

5. CONCLUSION 

We use a state-of-the-art cosmological reionization sim- 
ulation to quantify the topology of cosmological reion- 
ization, using genus curve measurements of the neutral 
hydrogen density distribution on > IMpc/i^^ scales. 
Overall, the reionization process begins in an inside-out 
fashion on scales, where higher density regions become 
ionized earlier than lower density regions. Reionization 
finally completes when the isolated neutral structures 
which remain are reionized in from outside-in. 

The sequence of topological transitions can be mapped 
to four phases: 

(1) pre-reionization: z > 15, (2) pre-overlap: 10 < z < 
15, (3) overlap: 8 < z < 10, (4) post-overlap: 6 < z < 8. 

In Phase (1) the genus curve is consistent with a 
Gaussian density distribution and topology remains un- 
changed as the universe evolves via linear evolution. In 
Phase (2) the number of HII bubbles increases gradu- 
ally with time, coincident with the increase in the am- 
plitude of the genus curve. Most of this occurs in the 
highest density regions which begin star-formation ear- 
lier. At the later stages of this phase, HII bubbles start 
to merge/percolate, resulting in a turn-over and sub- 
sequent decrease in the amplitude of the genus curve. 
An isolated Hll-bubble dominated topology is evident 
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in the genus curve during this phase. In Phase (3) the 
HII bubbles rapidly merge, manifested by the precipitous 
drop in the amphtude of the genus curve. The reioniza- 
tion of the IGM occurs rapidly as the ionizing photons 
are now free to stream freely between connected HII re- 
gions, increasing local ionizing fluxes. In Phase (4) the 
IGM is thoroughly ionized and the genus curve is consis- 
tent with an increasingly diminishing number of isolated 
neutral islands, with the rcionization of the universe at 

1 Mpc/i~^ complete by a redshift of ^ « 6. 

These results represent a concrete identification of the 
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